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ABSTRACT 

We compute the nonlinear development of the instabilities in C-shocks 
first described by Wardle, using a version of the ZEUS code modified to 
include a semi-implicit treatment of ambipolar diffusion. We find that, in 
three dimensions, thin sheets parallel to the shock velocity and perpendicular 
to the magnetic field lines form. High resolution, two-dimensional models 
show that the sheets are confined by the Brandenburg & Zweibel ambipolar 
diffusion singularity, forcing them to numerically unresolvable thinness. Hot and 
cold regions form around these filaments, disrupting the uniform temperature 
structure characteristic of a steady-state C-shock. This filamentary region 
steadily grows as the shock progresses. We compare steady-state to unstable 
C-shocks, showing excitation diagrams, line ratios, and line profiles for molecular 
hydrogen lines visible in the K-band, with the Infrared Space Observatory, and 
with NICMOS on the Hubble Space Telescope. 
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1. Introduction 

Within molecular clouds, interstellar gas becomes dense enough to recombine almost 
completely, reaching levels of fractional ionization less than \ = 10 -4 . Positively charged 
ions or grains carry the magnetic field in most regions except the central portions of 
protostellar disks, where neutral number densities exceed 10 10 cm -3 and electron-ion drift 
becomes important (e.g. Konigl 1989) . The charged particles carrying the field couple 
in turn with the neutral gas through collisions. When the fractional ionization is high, 
neutrals collide frequently with ions, and the field is coupled tightly to the neutrals. As the 
number of ions drops, they begin to resemble a sieve, and it becomes easier and easier for 
neutral gas to slip between adjacent ions, so that the field decouples from the neutral gas, 
diffusing through it. This process is called ambipolar diffusion, and was first understood to 
be astrophysically important by Mestel and Spitzer (1956) and Mouschovias (1979). In 
protostellar cores, ambipolar diffusion of the field out past infalling gas allows mass to fall 
in to the central region without dragging field along with it. In protostellar disks, ambipolar 
diffusion may allow the magnetic field to drive jets, reducing the angular momentum of 
accreting gas. 

One of the most directly observable effects of ambipolar diffusion is to spread out 
shock waves in dense gas (Draine 1980). Normally, the thickness of shocks is determined 
by gas dynamical or MHD dissipation in the shock front. Typical jump or J-shocks in 
the interstellar medium have thicknesses much shorter than the cooling length of the gas, 
so immediately behind the shock front, the gas gets quite hot, reaching the post-shock 
temperature predicted by the shock jump conditions for adiabatic gas. On the other hand, 
if a shock wave moves faster than the sound speed, c s , but slower than the Alfven speed in 
the rarefied ions, vm = |B| 2 / y/Aixpi, where B is the magnetic field and pi is the ion density, 
a compression wave can travel through the field and ions, moving faster than the neutral 
shock front. The accelerated ions can then accelerate the neutrals through collisional drag, 
spreading the shock front into a broad continuous shock, or C-shock, in which ambipolar 
diffusion maintains the thick shock front. C-shocks typically have thicknesses greater 
than the cooling length, so the gas cools as it is accelerated, and never reaches the high 
temperatures predicted for J-shocks. A recent review covering these points was written by 
Draine & McKee (1993). 

C-shocks were shown to be unstable by Wardle (1990, 1991a,b) . The mechanism 
depends on the field lines in a steady-state C-shock being perturbed so as to no longer be 
exactly perpendicular to the flow, as shown in Figure [I]. The ions have low inertia and 
collide frequently with neutrals, so a drag force due to the component of the neutral velocity 
parallel to the perturbed field lines drives the ions along the field lines into the valleys and 
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away from the peaks of the field lines. 

The density of ions in the valleys depends inversely on the wavelength of the 
perturbation for constant perturbation size, for two reasons. First, the ions have longer to 
travel to reach a valley at longer wavelengths, and second, the field lines are more normal 
to the flow at longer wavelengths, so the component of the neutral drag force driving ions 
into the valleys is lower. As the ions concentrate in the valleys, and rarefy at the peaks, 
the drag of the neutrals on the ions changes proportional to the ion density, transferring 
neutral momentum to the field, and producing a growing perturbation that leads in turn to 
more ion accumulation. 

There are stabilizing forces at both long and short wavelengths that act against the 
perturbation, due to magnetic pressure and tension, respectively. At long wavelengths, 
magnetic pressure acts on field lines compressed into valleys by the increased drag forces 
there, slowing or preventing further growth. At short wavelengths, magnetic tension acts 
against the curvature of field lines in valleys and peaks caused by increasing perturbation 
size, again slowing or preventing further growth. As a result, there exists an intermediate 
wavelength of maximum growth, computed by Wardle (1990) in his linear analysis. In 
three dimensions, the fastest growing modes were found by Wardle to be slabs, with no 
variation in the third direction. If there were variation, the compressed magnetic fields in 
the valleys could twist out sideways into the third dimension, releasing magnetic pressure, 
but decreasing the drag forces acting along the field lines, and so slowing the growth of the 
instability. 

Prior computations including ambipolar diffusion have focused on the collapse of 
protostellar cores and on the stability of C-shocks. The first dynamical computation to 
include ambipolar diffusion was done by Black & Scott (1982) using a two-dimensional, 
first-order code that included an iterative approximation to an implicit method of computing 
the ambipolar diffusion. Paleologou & Mouschovias (1983) described a one- dimensional 
code that included an adaptive mesh and an implicit method for the ambipolar diffusion. 
Fiedler & Mouschovias (1992) described a two-dimensional code with adaptive mesh and an 
implicit method. Mouschovias & Morton (1991) computed protostellar collapse problems 
by integrating over the vertical structure of an accretion disk and using an implicit method 
to model the horizontal, axisymmetric flow. Toth (1994, 1995) used a two-dimensional, 
flux corrected transport code to investigate C-shock stability, using both explicit and 
implicit methods. Mac Low et al. (1995) added ambipolar diffusion to the ZEUS codeQ 
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(Stone & Norman 1992a, b) and investigated its effects on Balbus-Hawley instabilities as 
a first application. 

In this paper we examine the nonlinear evolution of the Wardle instability with 
high-resolution two-dimensional and three-dimensional computations, extending the work 
of Mac Low et al. (1995) by adding the ion mass conservation equation, and a semi-implicit 
treatment of the ion momentum equation, following Toth (1995), in order to compute two- 
and three-dimensional models of unstable C-shocks. A preliminary version of this paper 
was presented by Mac Low & Smith (1997) . Stone (1997) describes similar computations. 

2. Scales and Parameters 

The time and length scales for our problem are the time for material to flow through 
a C-shock (e.g., Toth 1995; Wardle 1990), given by the time required for every neutral 
to collide with an ion tfi ow = t n i = l/jpnj, where the variables are defined below, under 
equation (f|); and the thickness of the C-shock, given by L S hk = \/2vAntfiow 

The time- dependent shock flow pattern is fully determined by four control parameters: 
the Alfven number A n , the Mach number M, the ion fraction and the initial field orientation. 
The Wardle instability is insensitive to the ion fraction and Mach number so long as they 
allow for the existence of a C-shock at all. Therefore only one parameter, A n , is necessary 
to describe the whole class of transverse, low-ionisation, cold flows: a single flow simulation 
with a fixed Alfven number is relevant to a wide range of conditions. 

The time scales for ionization and recombination can be important as well. Specifically, 
if the recombination timescale is significantly shorter than the flow timescale, the possibility 
exists of the instability being suppressed. If the ions are primarily molecular ions such 
as H 3 + , dissociative recombination rates may be high enough for this to be relevant 
(Flower & Pineau de Forets, personal communication). We will only consider this case in 
the extreme limit of constant ion density, which we will show does indeed suppress the 
instability. 

In this paper we discuss both two- and three-dimensional models. Our "standard 
model" has neutral Alfven number A n = (v s / B ) ^/Anpno = 25, initial number density of 
H atoms n# = 10 5 cm -3 , (note that p n = 1.4n#), fractional ionization x — 10~ 4 , shock 
velocity v s = 5 km s _1 , sound speed c s = 0.01 km s -1 , ion mass = 10m#, and neutral 
mass m n = 7/3w,h, assuming 10% He. We use an ion-neutral collisional coupling constant 
7 =< cxwq > I (mi + m n ) = 9.21 x 10 13 , where we take o = 10 -15 cm 2 , and an effective 
velocity of Wq = 1.9 x 10 6 cm s _1 . We quote our results in terms of the length, time, and 
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velocity scales given above, however, as this run will give results appropriate for A n = 25 
shocks in a wide variety of conditions. 

The resolution of our runs can best be expressed by giving the number of zones in the 
unperturbed shock thickness, L s hk- As shorthand, we use the format Lnn, where nn are 
the number of zones. We have performed runs with resolutions ranging from L60 to L240 
(note that Stone [1997] used a standard resolution of L21, though he focussed on slower 



shocks that had longer wavelengths of maximum growth, as we discuss in § [5.1| ). For our 
two-dimensional cases we used grids ranging from 160 x 40 to 640 x 160 zones, while for the 
three-dimensional cases we used grids with 160 x 40 x 40 and 267 x 80 x 40 zones. 



3. Numerical Methods 

Mac Low et al. (1995) described the basic interface with the ZEUS code. Summarizing, 
that work made four approximations: isothermality of ions, electrons and neutrals, no 
electron-ion drift, ion density dependent in power-law fashion on neutral density, and no 
ion inertia or pressure. This allowed us to neglect, respectively, the energy equation, Ohmic 
diffusion, and the equations of ion mass and momentum conservation. This approach has 
proved adequate for modelling protostellar disks, but broke down for both physical and 
numerical reasons when applied to C-shocks. The Wardle instability depends on the flow of 
ions along buckling field lines in the shock front, and so was suppressed by the neglect of 
the ion mass conservation equation. 

Neglect of ion inertia and pressure had allowed the replacement of the ion momentum 
conservation equation by an algebraic equation expressing the balance between Lorentz 
forces and ion-neutral drag in determining the drift velocity between ions and neutrals. This 
approach is physically accurate and allows time steps determined by the equivalent of the 
Courant condition for ambipolar diffusion, At < ir^fPip n (Ax) 2 /|B| 2 (Mac Low et al. 1995). 
However, both Toth (1994) and we have found that this approach is numerically unstable 
in the presence of steep velocity gradients as in C-shocks. Toth (1994) found that even 
an explicit treatment of the ion momentum equation was insufficient to entirely suppress 
numerical instability; Toth (1995) showed that a semi-implicit treatment of this equation 
was far more stable, an approach we have followed with some modifications described below. 

The equations we solve in the current version of the code are then the neutral and ion 
continuity equations, the induction equation, and the neutral and ion momentum equations: 



dp n /dt = -V-(p„v„) 
dpi/dt = -V • (pm) 



(1) 
(2) 
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p n (dv n /dt) = -p n (v n • V)v n - VP„ 

+7PiPn(Vj - v n ) (3) 

Pi(dvi/dt) = -pi(vi- V)vj- VPi + 7Pip n (v n -Vi) 

+ (1/4tt)(V x B) x B (4) 

(dB/dt) = Vx(v.xB) (5) 

V-B = (6) 

where the subscripts i and n refer to the ions and neutrals, p, v, and P are density, velocity 
and pressure for each fluid, B is the magnetic field, and 7 is the collisional coupling constant 
between the ions and neutrals. 

We treat the ions as a separate fluid in the code, using the standard ZEUS algorithms 
to update them, except in the momentum equation where we make several changes to 
implement the semi-implicit algorithm. The basic idea of this new algorithm is to treat 
the stiffest terms in the equation implicitly, but maintain the standard operator splitting 
for the other terms. Toth (1995) included the drag term and the transported fluxes in his 
semi-implicit computation. We find that we also needed to include the magnetic pressure 
forces in order to maintain balance between them and the drag terms. 

We rewrite equations @ and (|J) in conservative form using momenta s = pv, and 
momentum fluxes F n and Fj, and place them in finite difference form. We then get two 
equations in the two unknowns s™ +1 and s™ +1 for our implicit system: 



s n+lyn+l _ ^nyn 
At 



= -VF n (7) 

+ 7 ^ n+1 (Pnsr fi -p l sr i ), 



s n+lyn+l _ s nyn B 2 



At v 8tt 



where time levels are given by superscripts, the finite difference time step is At, and the 
inclusion of zone volumes V allows for a moving grid. We find the semi-implicit update to 
the momenta by solving these equations for s™ +1 and s" : 

s^ +1 = (9) 
|(s™K"-VF n At)(l + 7 Atp„) 



+ 7 Atp n (s^ n -V 

{[l+7At(p n + p J )]K" +1 }- i 



Fi+ — 

07T 



At) 
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sr +1 = (io) 

+ jAt Pi {slV n - VF n At) } 
{[l + 7 At(p n + p J )]^ +1 }- 1 . 

Of course, the magnetic pressure updates must be removed from the source step where they 
would normally be performed. 

Unfortunately, including the ion momentum equation limits our timestep to 
At < Ax j v Aii and we can no longer use the substepping scheme described by Mac Low et 
al. (1995). As a result, computations through a flow time require many xlO 5 cycles. The 
total number of cycles goes as x 1 ^ 2 -, however, so increasing the ionization fraction helps, 
so long as it does not get so high that no C-shock is possible. Toth (1995) showed that a 
faster algorithm could be implemented by going to a fully implicit method, at the cost of 
greatly increased code complexity. 

For initial conditions we used a numerical solution for an isothermal C-shock (Smith 
& Mac Low 1997) interpolated onto our grid. We then added a random zone-to-zone 
fluctuation to the neutral velocity. In order to maintain the shock stationary on our grid, 
we set up the grid with defined values for the field variables at both ends — that is with 
"inflow" boundaries, but with the right side having a negative inflow velocity. We used 
periodic boundaries along the sides of the grid. 

We can maintain the steady-state C-shock indefinitely in one dimension (Smith & Mac 
Low 1997). In two and three dimensions, we can maintain it for many thousands of cycles 
until the physical instability sets in, as shown in Figure [| We also showed that we converge 
on the steady-state solution starting from a discontinuous initial condition in Smith & Mac 
Low (1997). 



F i + — 
87T 



A*)(l+ 7 Atpi) 



4. Instability 

4-0.1. Linear growth 

Let us begin by verifying that we are indeed seeing the instability in the linear regime. 
We follow Stone (1997) in using the magnetic energy in the field component parallel to the 
shock velocity, B%, as a diagnostic. The growth rate depends almost entirely on the neutral 
Alfven number of the gas and the angle of the shock to the field, as discussed by Wardle 
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(1990) and Toth (1994). The analytic growth rate s can be derived from Figure 7 of Wardle 
(1990). Note, however, that that Figure gives the quantity stfi ow , and that it appears to 
assume the value of tji ow given in Wardle's equation (2.12). The parameters quoted by 
Wardle would actually give a coefficient of 5.4 rather than 5.0, but that does not agree with 
our numerics, while 5.0 does. (Note that, for the parameters we use, the coefficient is 6.5, 
but that is due to our different choices of and so forth.) 

In Figure ^| we show the parallel field energy as a function of time for our standard 
case for numerical resolutions of L60, L120, and L240. For this case, we derive a value 
of s = 42.5£/; OI1) ~ 1 , using Wardle's value of the coefficient. In Figure |3] we show curves 
corresponding to s — (41, 42.5, 44)i// ou) -1 , showing the excellent agreement with the analytic 
value. We have also measured the growth rate for our A n = 10 run, and get equally good 
agreement with the analytic value of s — A.2tfi ow ~ l . 

4-0.2. Nonlinear Development 

The striking feature of the nonlinear development is the formation of long, thin 
filaments of dense gas. These filaments begin to merge with each other so that the 
wavelength of the instability grows as it becomes increasingly nonlinear. (The nonlinear 
Rayleigh- Taylor instability also develops merging filaments, as described by Read [1984] 
and Youngs [1984] ). We show this development in Figure |j. 

In our standard model, numerical instabilities prevent computation of the further 
development of the filaments after the last time shown in Figure However, we have been 
able to extend our computations of lower neutral Alfven number shocks to later times, as 
shown in Figure |5[ We confirm the finding of Stone (1997) that, in the nonlinear stage 
when the distance between filaments reaches the shock thickness, a new filament appears 
between the existing ones. 

The thickness of the region between the initial ion deceleration zone and the final 
post-shock region increases with time as shown in Figure |6|. The thickness grows with 
a velocity just equal to the post-shock velocity of the shock. It appears that once the 
instability is established, it does not relax back to a uniform medium easily, if at all. The 
post-shock state is thus a tree of merging filaments (or in three dimensions, merging sheets, 
as we discuss below.) 

We find that the filaments' width drops to the single zone level at every resolution that 
we have run the problem. Stone (1997) suggests that gas pressure eventually limits their 
thickness, but points out that the resolution needed to see this effect would be enormous, 
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requiring a dynamic range over 10 4 . We have confirmed this empirically by running several 
different Mach 25 runs with sound speeds as high as OAvAn at our highest resolution, 
with no discernible difference in their evolution. However, as we show next, the filaments 
are probably manifestations of the Brandenburg & Zweibel (1994, 1995; hereafter BZ) 
ambipolar diffusion singularity, which they showed in the simplest case overwhelms thermal 
pressure, and so may well become extraordinarily thin. 

4.1. Ambipolar Diffusion Singularity 

In Figure [7] we show the components of the magnetic field parallel and perpendicular 
to the shock velocity, along with the magnetic pressure. Although the perpendicular 
component remains dominant, as can be seen in the magnetic pressure plot, the parallel 
component undergoes a sharp reversal at each filament. A profile of the parallel field is 
shown in Figure [8[ The field lines thus are basically perpendicular, but with extremely 
sharp kinks in them at each filament. The reversal in the parallel component appears to be 
enough to drive the BZ ambipolar diffusion singularity. This singularity steepens magnetic 
pressure gradients in the presence of ambipolar diffusion, due to the nonlinear nature of the 
ambipolar diffusion. 

The actual field behavior that we compute can be understood analytically. The 
simulations demonstrate that the magnetic field parallel to the flow, B x , grows exponentially 
(Fig. |3p, evolving towards a saw-tooth structure as the filaments become thinner. That 
is, B x becomes a linear function of y. In contrast, the ambipolar diffusion analysis of 
Brandenburg and Zweibel (1994) delivers a relation B x oc y 1 ^ 3 . Their analysis, however, 
assumed the magnetic field to be zero within the filament whereas here we have the 
opposite case: the field B y dominates. Hence to understand these different relationships we 
undertake a more general analysis. 

The aim, following BZ, is to search for hydrostatic steady flow solutions to the 
induction equation. To begin we take the streaming (or drift) velocity from equation (|J) 

v fl = (v„ - v<) = -(l/47r 7 p^ n )(V X B) X B (11) 

to eliminate the ion velocity from the induction equation (|5]): 

(dB/dt) = V x (v n x B) + (l/47r7pip n )V x [{(V x B) x B} x B] . (12) 

We take a uniform flow in the x-direction i.e. B x (y), v xn (y) and B y = Bq, a constant as 
required for vanishing divergence. Ignoring field amplification due to the first shear term 
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(that is, taking v xn constant), the induction equation becomes 



dB x 

'~dT 



d 2 B x 
dy 2 



+ 2B X 



' dB x 
, dy 



(13) 



The solution to the steady state form is 

y = k 1 + k 2 



n tan 3 9 

2 tan 9 H T - 

sin"' 9 



(14) 



where the field angle 9 = arctan (B x /B ) , and k\ and ki are integration constants. The 
second term dominates and yields the Brandenburg & Zweibel (1994) solution in the limit 
B — > 0. Such a solution may be relevant to the study of parallel-field C-shocks that we 
plan. The first term, however, dominates for the case studied here, and indeed predicts a 
linear y 



B x behaviour. 



This singularity steepens magnetic pressure gradients in the presence of ambipolar 
diffusion, due to the nonlinear nature of the ambipolar diffusion. Among other things, the 
resulting steep gradients trigger numerical instabilities, visible as the single-zone regions of 
opposed field along the filaments. These regions produce anomalous fields and extremely 
high ion velocities, reducing the timestep to impractically small values. Physically, the 
exceedingly sharp discontinuities in field direction produced by this singularity will drive 
strong current sheets along the filaments, with uncertain consequences probably including 
heating and possibly the production of energetic particles. 



4.2. Three-Dimensional Structure 

To study the three-dimensional structure of the instability we performed two runs 
at L120. One was of our standard case, while the other had a zone-to-zone perturbation 
5v n i/v n x = 10~ 6 rather than 0.01. We found that the thin filaments seen in the two- 
dimensional models developed in the low perturbation run into coherent mildly rippled 
layers, as shown in Figure ^. In the standard case, on the other hand, we found that each 
vertical plane developed essentially independently, resulting in spikes rather than layers. 
This appears to have been an effect of the initial perturbation, suggesting that the structure 
in the third dimension in nature will be determined by the existing perturbation structure 
in the pre-shock gas unless it is quite uniform on scales of L^. 
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4.3. Resolution Studies 

In order to understand the effects of numerical diffusivity on our results, we ran 
our standard case in two dimensions at resolutions of L60, L120, and L240, and in three 
dimensions at resolutions of L60 and L120, as we discussed at the beginning of this section. 
We show the density distribution at the end of the run for these three different cases in 
two dimensions in Figure [10|. Examining this figure, it is clear that the filamentary nature 
of the nonlinear instability has been well resolved, but that we have not converged on the 
exact value of the wavelength of maximum growth. As was shown in Figure |6], the thickness 
of the shocked layer is converged to within a few percent. 

The thickness of individual filaments in the transverse direction appears to be 
determined entirely by the dissipation scale physically, due to the BZ ambipolar diffusion 
singularity, as we discussed above, in § |4.0.2| Effectively, they become infinitely thin. Indeed, 



they reach thicknesses of one or two zones at every resolution, at which point the extreme 
gradients cause numerical instabilities. 



5. Variations on the Theme 

5.1. Mach Number Dependence 

In the linear regime, the growth rate and wavelength of maximum growth depend on 
the neutral Alfven number, A n , as shown by Wardle (1990), and confirmed in some detail 
by Toth (1994, 1995), as well as in Figure ^. The slower growth rate and longer wavelength 
of maximum growth in the shock with A n = 10 shown in Figure as compared with 
our standard model, agree with the linear theory. We find qualitative differences in the 
nonlinear development of the instability at different A n as well, as can be seen by comparing 
Figures f| and [5[ At lower values of A n , the larger wavelength of maximum growth means 
that there is more space between the resulting filaments, and each one is thicker compared 
to the initial shock thickness. 



5.2. Ionization Equilibrium 

Until now, we have assumed that ionization and recombination timescales were long 
compared to the flow time through the shock tfi ow . If both ionization and recombination 
timescales are very short compared to the flow time, the ion density will remain constant 
with time, or vary as a power law of the neutral density, suppressing the instability. The 
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version of ZEUS described by Mac Low et al. (1995) works in this regime, so we ran our 
standard case with constant ion density using that version. We did indeed find that after 
0Atfi ow there was no indication of the Wardle instability. At that point, however, the run 
was ended by numerical instabilities caused by the large ion-neutral streaming velocities, 
which caused strong zone-to-zone fluctuations in the computed parallel velocities. In 
Figure |TT] we compare the neutral density distribution for our standard initial conditions to 
the result for constant ion density, to show the lack of evolution. 



5.3. Initial Perturbations 

The form and strength of the initial perturbations we use to trigger the instability can 
make a difference in the results. In our first runs we used a 1% zone-to-zone perturbation 
in density; later we followed Toth (1995) in using perturbations in parallel neutral velocity 
v n i instead, with magnitudes of 1% and 0.0001%. We found that 1% perturbations were 
already large enough to significantly change the final nonlinear results, especially in three 
dimensions, as we discussed in § [4.2| . Even in two dimensions, the larger perturbations 
change the wavelength of maximum growth significantly, as shown in Figure O, where the 
nonlinear development of the standard case is compared to that of a low perturbation case, 
which develops a longer wavelength. 



6. Implications for Molecular Hydrogen Emission 

One measurable consequence of the C-shock instability is that the line emission directly 
from the shocked layer will differ from the steady-state case, and should even be a slowly 
varying function of time. Warm C-shocks vibrationally excite the hydrogen molecules. We 
may thus expect that the fraction of H2 in higher vibrational levels would decrease with 
time as frictional heating is spread out over a thicker layer. To test this, we employ the 
"cool C-shock" approximation, that requires the temperature T <C 2rriHV 2 s /k. This is valid 
provided cooling is very efficient so that the thermal pressure and pressure gradient can be 
neglected. Then, the isothermal assumption is not critical to the flow pattern so long as a 
low temperature is chosen. In this case, an effective temperature can be derived at each 
location by equating the frictional heating to the molecular and atomic cooling (see Smith 
& Brand 1990a, Smith 1993). 

The cooling functions we use here have the form n^'T 01 erg s _1 cm -3 . We study 
two cases: the high-density regime with p n ~ 10 8 cm -3 , where H2O cooling dominates, 
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/ = 1, and a = 1.2; and the low density regime with p n ~ 10 4-6 cm -2 , where H 2 cooling 
dominates, / = 0, and a = 3.3 (Smith 1993). Frictional heating is oc n n riivf n . This yields 
T a = t]vf n ni/nl, where we choose the constant rj so that the initial C-shock possesses 
a maximum temperature of 2,500 K. The temperature distribution for the two different 



cooling cases is shown in Figure [13] for the initial state and a time 0.36tfi ow , after the 
instability has gone nonlinear. 

The column density of H2 molecules in an excited state j with energy level kTj is given 

by 

= \ gj N(R 2 ) exp(-r i /r)/Z(T). (15) 

The statistical weights gj = 4ip(2Jj + 1), where if) is the fraction of H 2 molecules in the 
ortho state (assumed to be 0.75 here), Jj is the rotational quantum number, and Z(T) is 
the partition function 

Z{T) = 6.13510" 3 T [1 — exp(-5850/T)r 1 , (16) 

(Smith, Davis & Lioure 1997), assuming local thermodynamic equilibrium. The intensity of 
a line is then given by 

Lj = (hc/X^NjAj, (17) 

where the radiative deexcitation rates Aj (Turner, Kirby-Docken, & Dalgarno 1977), along 
with statistical weights, wavelengths (Black & Dishoeck 1987), and excitation temperatures 
for a number of important lines are given in Table 1. 

Neufeld & Stone (1997) concluded that the Wardle instability would not markedly 
change the emission expected from steady-state C-shocks. However, they only considered 
shocks with neutral Alfven number A n = 10 and n Q = 10 5 cm -3 , which will be dominated 
by molecular hydrogen cooling. While we agree that the emission from these shocks is not 
strongly modified by the Wardle instability, we now show that shocks in other regions of 
parameter space behave much differently. 

One common method of deducing molecular excitation is examining the ratio of the 
strong H 2 lines 1-0 S(l) and 2-1 S(l), which should increase as a gas cools and the second 
vibrational level becomes less populated. We calculate the column density of molecules in 
each grid zone, and then sum over zones to find the total number of molecules on the grid in 
each upper energy level, Nj. Radiative transfer could be added in this treatment. However 



these quadrupole transitions are expected to be optically thin. Figure [L4] shows the ratio of 
line intensities L(l-0 S(l))/L(2-1 S(l)), ignoring foreground extinction, as a function of time 
and numerical resolution for both our standard case and our model with A n = 10. We find 
that the instability can shift the ratio strongly when H 2 cooling dominates at A n = 25, 
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and more weakly in other cases. We have not converged on the actual value of the shift 



at our current numerical resolution, primarily because the hot filaments seen in Figure [13 
are not yet well resolved. As unincluded physics is probably also important there, higher 
resolution runs are not particularly justified. 

A more general method of determining the molecular excitation is to correlate column 



densities (which can be derived from line intensities using equation [17| and Table 1) 



in a number of lines against the excitation temperatures of the lines. Because of the 



strong temperature dependence of the column densities in equation [L5|, we scale the 
column densities by gj exp(— Tj/2000K), the temperature dependence of a slab of gas with 
temperature 2000 K, where the statistical weights gj and excitation temperatures Tj for 



important lines are given in Table 1. In Figure ITH we show the resulting time-dependent 



development of the excitation conditions for our standard case with either H2O and H2 
cooling dominant. The instability forms both hot and cold regions compared to the 
steady-state C-shock, and so produces a marked shift in the integrated excitation conditions 
across the shock. 

Finally we can examine the detailed line profiles produced by stable and unstable 
C-shocks. We compute them by binning the grid into angle-dependent velocity bins, and 
then adding up the intensity of emission from all zones in each bin (scaled by the square 
of the zone-size for runs of different linear resolution) to get a spectrum. We smooth the 
resulting profiles with a 1 km s _1 Gaussian (except for the narrow profile from a line of sight 
perpendicular to the shock velocity, which we smooth with an 0.1 km s -1 Gaussian). The 
intensity scales displayed are arbitrary, but consistent across all plots to allow comparison 
between lines and angles. 

We first examine the effects of angle and numerical resolution on profiles of the bright 
1-0 S(l) line in Figure |TB|. We show the line profile for lines of sight parallel to the shock 
velocity, at 45° to it, and perpendicular to it, for shocks dominated by H 2 and H 2 cooling. 
Note that the velocity scale is different for each angle in the figure. We find that the line is 
broadened to the Alfven velocity even on lines of sight normal to the shock velocity. We find 
that the details of our profiles are dependent on our numerical resolution, as can be seen by 
comparing different line profiles in each plot, which are from runs differing by a factor of 
two in linear resolution. However, the gross shape and intensity of each line appears well 
described. The smaller scale variation is due to the dependence of the results on the details 
of the filament structure where most of the hot gas is produced. Again, the structure of 
these filaments is also liable to be influenced by physics not yet included in our model. 

Finally, we examine the profiles of a number of lines observable with ISO (Figure [17] ), 



in the K-band, or with NICMOS on the Hubble Space Telescope (Figure [18]), for cooling 
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dominated either by H 2 or by H 2 0. We show in grey the profiles predicted for a steady-state 
C-shock viewed parallel to the shock velocity. These are very well resolved, and represent 
the first publication of predicted line profiles for most of these lines; they agree well with 
the general predictions of Smith & Brand (1990b) . In black we show the same profiles 
at the end of our standard run, showing the shifts in line profile and intensity due to the 
growth of hot and cold regions due to the instability. We also give the derived total column 
densities from these lines in Table 1 to allow computation of arbitrary line ratios. The 
convergence of the unstable profiles can be judged by comparison with Figure [Tj| The 
general direction and size of the shift from the steady-state value is probably correct, but 
the quantitative results are probably good to no better than 20%. The general trend we 
observe is a strengthening of lines at the soft and hard ends of the spectrum due to the 
broader range of temperatures in the C-shock, as well as changes in the velocity structure 
due to the filaments. 



7. Summary 

• The fastest growing modes of the Wardle (1990) C-shock instability are layers with 
thickness of order 0.1 L s hk in the plane parallel to the field and the shock normal, 
in agreement with Wardle's linear analysis. As predicted, our standard shock with 
neutral Alfven number A n = 25 goes unstable quite rapidly, becoming very non-linear 
in a fraction of tfi ow . 

• When the shock goes non-linear, the layers collapse into thin, dense filaments (sheets 
in three dimensions) with density more than a factor of two above the expected 
post-shock density. A dense sheet also forms across the shock front at the downstream 
end of the filaments. These filaments are physically probabably too thin to be 
well-resolved numerically. Toth (1994, 1995) did not see these filaments, possibly due 
to the diffusivity of his code, but Stone (1997) finds a similar morphology. 

• As the instability becomes more nonlinear, the sheets collapse toward each other, 
as seen in the upper panels of Figure [|. Thus, longer wavelength modes become 
dominant as time goes on. This process has proceeded in some of our runs until only 
one wavelength is left on the grid. Rayleigh- Taylor instabilities undergo a similar 
lengthening of the dominant mode in the non-linear regime. New filaments appear 
when the distance between filaments becomes comparable to shock thickness. 

• The thickness of the shock front increases linearly as the filaments grow, with the ion 
velocity dropping sharply at the ends of the growing filaments. 
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• The Brandenburg & Zweibel (1994, 1995) ambipolar diffusion singularity steepens 
the fields along the filaments, forcing the filaments to be ever thinner, until the 
code breaks down due to the resulting sharp gradients. Strong current sheets form 
along the filaments as a result, with consequences that probably include heating and 
ionization, and possibly the production of energetic particles. The physics of these 
filaments have not been completely captured with our approximations. 

• We compare the integrated excitation conditions across stable and unstable C-shocks, 
giving comparisons of scaled column density to line excitation temperature, line 
profiles, and line ratios for our standard case. These results are sensitive to the as 
yet unmodeled physics of the filaments, as well as numerical resolution of the very 
thin filaments, but probably give a generally correct description of the effect of the 
instability on observables from C-shocks. 

• Significant variations in emission line properties are predicted on timescales of 
fractions of tfi ow = (30/nj) yr, where rii is the ion number density. Thus, multiple 
observations of suspected C-shocks in dense regions should prove rewarding. 



We thank J. Stone for generously sharing his results prior to submission, S. Beckwith 
for emphasizing the need for predictions of multiple line properties, and E. Zweibel for 
discussions of the ambipolar diffusion singularity. MDS thanks the DFG for financial 
support. Computations were performed at the MPG Rechenzentrum Garching. 
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Table 1. Representative Molecular Hydrogen Lines 
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Table 1 — Continued 
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Fig. 1. — The Wardle instability occurs when field lines are perturbed slightly. The neutral 
drag forces ions into the valleys, increasing the drag in the valleys and reducing it at the 
peaks. As a result, the valleys are further deepened, while magnetic pressure expands the 
peaks as well, driving an exponential instability. 

Fig. 2.— Ion (earlier drop) and neutral (later drop) velocities scaled to the shock 
velocity. Dotted lines show the analytic solution, while solid lines show the two-dimensional 
higher resolution (grey) and three-dimensional (black) solutions after 0.11tfl O w The three- 
dimensional solution has been shifted by 1.5L shk for visibility. 

Fig. 3. — Growth of instability, as indicated by the growth of magnetic energy in the field 
component parallel to the shock velocity, B^, normalized to the initial magnetic energy on 
the grid By . Our standard case in two dimensions with neutral Alfven number A n = 25 is 
shown at resolutions of L60 (dotted line), L120 (dashed line), and L240 (thick solid line). 
Overlaid are three exponentials (thin solid lines) of form e st , with s = (41, 42.5, 44)t// oll ,~ 1 , 
where the analytic value is s = 42.5tfi ow ~ 1 , as discussed in the text. In three dimensions 
we show a run with the same parameters at a resolution of L120 and with smaller initial 
perturbations (thick solid line). The overlaid exponential (thin solid line) in this case is the 
analytic value of the growth rate, normalized appropriately. 

Fig. 4. — Time history of (a) neutral parallel velocity, (b) neutral density, and (c) ion density 
for our standard model with a neutral Alfven number A n = 25, and other parameters given 
in the text, run at a resolution of L240 on a grid of 640 x 160 zones. 

Fig. 5. — Time history of ion density for our model with neutral Alfven number A n = 10, 
run at a resolution of L120 on a large grid of 400 x 240 zones. The first indications of a 
density enhancement where another filament will form can be seen in the last panel, between 
the upper two filaments. 

Fig. 6. — Shock thickness as a function of time at numerical resolutions of L60 (dotted), 
L120 (dashed), and L240 (solid), showing the near-linear growth in the size of the disturbed 
shocked region. 

Fig. 7. — Parallel and perpendicular components of the field at the end of our standard run. 
Note that the field remains dominated by the perpendicular component. 
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Fig. 8. — Profile across filaments of the parallel field component B x at the end of our standard 
run, at the midpoint of the grid in the x-direction. Values in individual grid zones are shown 
by triangles. Note the linear increase in the field, and the extremely sharp transitions in the 
filaments. 



Fig. 9. — (a) Comparison of neutral x-velocity distribution in the x — y plane for two- and 
three-dimensional versions of the standard model at a resolution of L s hk = 120 zones (total 
grid 267 x 80 x 40 zones), (b) Cuts through the y — z plane of the three-dimensional model 
at locations indicated by black lines in (a), scaled to the same palette. 



Fig. 10. — Comparison of results for runs with resolutions of L60, L120, and L240. The log 
of the neutral density distribution is shown for three runs with total grids of 160 x 40 zones, 
320 x 80 zones, and 640 x 160 zones, for our standard model. Note the convergence of the 
morphology and the wavelength of maximum growth, as well as the consistent increasing 
trend in the peak density, which has not yet converged. 



Fig. 11. — Comparison of runs with conserved ions versus constant ion density, showing that 
if ionization and recombination occur fast enough to maintain a constant ion density, the 
instability is suppressed. Log of neutral density is shown, with the same scaling across both 
images. The run with conserved ions is our standard model at a resolution of L shk = 120 
zones, while the run with constant ion density was run with ion fraction x = 10~ 6 for speed. 
(This should make no physical difference, aside from changing the absolute time and length 
scales.) The run at constant ion density had to be halted due to numerical instability at this 
time, but no sign of physical instability was seen. 



Fig. 12. — Density distribution for runs with the given initial perturbation strength (in 
both cases, perturbations consisted of zone-to-zone random changes in the parallel neutral 
velocity). The strong perturbation case is shown at a time t = 0.36tj/ ow , while the weak 
perturbation case is shown at a time t = 0.69t// ow . The strong perturbation standard case 
was only run with half the number of zones in the y-direction, so it has been doubled for 
comparison to the 640 x 320 zone weak perturbation case. 
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Fig. 13.— When drag heating balances radiative cooling, the equilibrium temperature 
can be computed assuming that either (a) water or (b) molecular hydrogen dominates the 
cooling. We show the log of temperature for our standard model before the instability sets in, 
and after it has become nonlinear at a time 0.36tfi ow . The earlier time shows the standard 
picture of C-shocks, while the later time emphasizes the difference the instability makes. 
Note the existence of warm gas at the temperatures predicted by the standard picture, but 
mixed with colder and hotter gas. 

Fig. 14. — Ratio of emission intensity from the 1-0 S(l) and 2-1 S(l) lines of molecular 
hydrogen for cooling dominated by (a) water and (b) molecular hydrogen for our A n — 25 
run. Lighter grey lines show a numerical resolution of L60, darker grey lines are L120, and 
black lines are L240. (c) The same ratio is shown for our A n = 10 run for both molecular 
hydrogen cooling (dashed) and water cooling (solid). 

Fig. 15. — Variation of excitation conditions in an unstable C-shock at times of (lightest 
grey), 0.19tfi ow , 0.29t// om , and 0.36t// ow (black) for cooling dominated by (a) water and (b) 
molecular hydrogen. The column density is scaled by gj exp(— Tj/2000 K), the functional 
form of column densities from a uniform slab of gas at a temperature of 2000 K, where the 
statistical weights gj and excitation temperatures Tj are given in Table 1 for many common 
lines. Thus, such a uniform slab of gas would produce a flat line in this plot. 

Fig. 16. — Variation with angle to the shock velocity of the line profile from the 1-0 S(l) line 
of H 2 at 2.121 (ttm in the final state of our standard case at resolutions of L60 (lighter grey), 
L120 (darker grey), and L240 (black), demonstrating the degree of numerical convergence 
expected for cooling dominated by (a) water and (b) molecular hydrogen. Note that the 
velocity axes are scaled by different values at each angle, with v s = 25v An for this shock with 
A n = 25. 

Fig. 17. — A comparison of line profiles from an unstable C-shock (black) to profiles from 
a steady-state C-shock with our standard initial parameters (grey) for a sample of H2 lines 
observable by ISO, assuming (a) water cooling and (b) molecular hydrogen cooling. The 
unstable C-shock profiles are computed from the final state of our highest resolution L240 
run. 
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Fig. 18. — A comparison of line profiles from an unstable C-shock (black) to profiles from a 
steady-state C-shock with our standard initial parameters (grey) for a sample of near-IR H2 
lines observable in the K-band, or with NICMOS on the HST, assuming (a) water cooling 
and (b) molecular hydrogen cooling. The unstable C-shock profiles are computed from the 
final state of our highest resolution L240 run. 



Wardle Instability 





0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 

x 1 L shk 



This figure "fig4a.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "fig4b.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "fig4c.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "fig5.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "fig7.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



2 




0.0 0.2 6% shk 0.6 0.8 




0.0 0.2 (ffif shk 0.6 0.8 




0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 
Z/L shk Z/L shk Z/L shk Z/L shk Z/L shk 



This figure "figl0.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



y / L shk y / L shk 

oooooo oooooo 

— ^ N3 CM 4^ Cn CD — ^ N3 CnI 4^ Cn 01 

OOOOOO oooooo 




I 



This figure "figl2.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "figl3a.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



This figure "figl3b.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9703172vl 



12 




12 




600 



400 



1 

■v. 

a 



200 



0-0 S(1) 
17.035 ftm 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



80 

2 60 
1 

■e 

£ 40 

c 

u 

c 20 



0-0 S(2) 
12.279 /urn 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



15 



i 

3 
o 

o 



10 



/\ 




o-o s(4) : 




8.025 fxm - 







-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



C 

o 
e 



20 
15 
10 



0-0 S(5) 
6.910 /um 




5 


-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



2.5 

I" 5 

•s 1.0 
a 

S 0.5 
0.0 



0-0 S(6) 
6.109 fim 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



50 
| 30 



•S 20 
c 
m 

3 10 




0-0 S(7) 
5.511 flm 



1.0 



9 0.8 



t 0.4 
c 
u 

S 0.2 
0.0 




0-0 S(8) 
5.053 fxm 





1.2 




1.0 








0.8 


I 






0.6 


>■, 




9 


0.4 


_£ 






0.2 




0.0 




0-0 S(9) 
4.695 fj,m 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



1500 



1 1000 

t 




0-0 S(1 
17.035 fjur 








1 500 





I J 



200 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



-P 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



50 

1 30 

>, 

■| 20 1 
d 

3 10 




0-0 S(4) 
8.025 /xm 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

V / v„ 



50 

9 40 
x, 30 

is 

■« 20 



•9 10 



0-0 S(5) 
6.910 yu,m 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

V / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 






1.2 : 




1.0 






"5. 
3 


0.8 - 


1 






o.6; 






tens 


0.4 : 







0.2 
0.0 



0-0 S(8) 
5.053 yum 






1.4 




1.2 


1 


1.0 






1 


0.8 






>, 


0.6 


a 






0.4 








0.2 




0.0 









0-0 S(9): 




i/\ 4.695 yum 







-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 





-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 




0.000 



-0.2 0.0 



0.2 0.4 

v / v 



0.6 0.8 





0.12 : 




0.10; 






1 


0.08 ■ 






i 


0.06; 






c 


0.04; 








0.02- 




0.00 




2-0 S(5) 
1 .085 fim 



-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



1.0 



0.40 



2 0.30 



0.20 



0.10 



0.00 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.040 



0.030 



0.020 



0.010 



0.000 




0.15 



-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 







1 


0.10 






i 








tens: 


0.05 


a 





0.00 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



2 3 





-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.060 
0.050 




0.000 



-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.050 



- 0.040 




0.10 



0.000 



-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.20 



2 0.15 



0.10 



0.05 



0.00 




-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.015 






0.010 




/ \ 3-2 S(3) 
/ y.201 fj.m 


0.005 






0.000 







-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



0.025 



0.020 



0.015 



■f 0.010 
c 

" 0.005 



0.000 





/ \4-3 S(3) _ 




/ C.344 fj,m 


/ 





-0.2 0.0 0.2 0.4 0.6 0.8 

v / v„ 



1.0 



